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ABSTRACT 


A procedure  is  described  to  determine  lateral  buckling 
loads  for  initially  straight  beams.  Loading  and  beam  geom- 
etry may  vary  along  the  length.  Warping  rigidity  is  not 
considered.  End  conditions  of  considerable  generality  may 
be  treated.  The  algorithm  depends  on  solving  a convergent 
sequence  of  sixteenth  order  eigenproblems . A computer  pro- 
gram implementing  this  procedure  has  been  developed  and  is 
described  herein. 
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The  problem  of  analysis  logically  divides  itself  into 
three  cases. 

1. )  The  case  where  warping  rigidity  may  be  neglected. 

2. )  The  case  where  warping  rigidity  is  nonzero 

throughout  the  entire  beam. 

3. )  The  case  where  warping  rigidity  is  zero  over  a part 

of  the  beam  and  nonzero  over  the  remainder. 

As  will  become  clear,  these  three  cases  must  be  treated 
separately.  The  third  is  most  difficult  and  has  not  been 
treated  in  the  literature.  The  second  has  been  addressed  by 
J.  E.  Brock;  cf.  reference  (1).  The  first  case,  where  the 
warping  rigidity  is  zero  or  is  neglected,  is  not  a special 
instance  of  the  second  and  a separate  analysis  is  required. 
Such  an  analysis  forms  the  subject  of  this  thesis.  The  de- 
velopment is  parallel  to  that  of  Brock  Cl)  to  an  extensive 
degree  and  many  of  the  figures  and  equations  herein  are  iden- 
tical to  those  in  reference  (1). 


II.  PROBLEM  SPECIFICATIONS 


The  coordinate  system  used  is  a right-handed  orthogonal 
x,y,z  system  with  positive  senses  as  indicated  in  figure  1. 
The  longitudinal  centroidal  axis  of  the  beam  lies  along  the 
+z  direction.  The  principal  axes  of  a normal  cross  section 
are  originally  in  the  x and  y directions  but  after  torsional 
displacement  (rotation)  they  are  in  the  £ and  n directions, 
cf.  figure  2 (page  10). 


V 


Figure  1.  General  configuration  and  coordinate  system. 

The  beam  displacements  are  defined  by  three  variables: 
vertical  deflection  v,  horizontal  deflection  u,  and  angular 
rotation  Positive  senses  for  v and  u are  as  indicated  in 
figure  1.  The  positive  sense  for  <fi  is  indicated  in  figure  2. 
The  deformations  u,  v and  <p  are  considered  to  be  "small"  in 
a mathematical  sense,  so  that,  for  example  sin  <t>  « <<>. 


Figure  2.  Positive  sense  of  0 and  source  of  bending 
moment  N-<f>M. 

Defining  the  problem  requires  specifications  for  five 
physical  characteristics  of  the  beam.  They  are  beam  length 

v 

L,  flexural  rigidity  for  bending  in  the  vertical  plane  Ell, 
flexural  rigidity  for  bending  in  the  horizontal  plane  EI2  , 
torsional  rigidity  C,  and  warping  rigidity  Cl.  The  last 
four  may  vary  along  the  length  of  the  beam. 

Torsional  rigidity  C is  defined  by  the  following  equa- 
tion. 

C = (TORQUE) /(ANGLE  OF  TWIST  PER  UNIT  LENGTH)  (1) 

ft 

For  simplicity  in  notation,  we  use  numerical  suffixes 
rather  than  numerical  subscripts . 


Warping  rigidity  Cl  is  related  to  nonuniform  torsion.  A 
discussion  here  of  this  property  would  be  too  digressive; 
cf.  section  5.3  of  reference  (2). 

It  is  the  underlying  purpose  of  this  entire  analysis  to 
determine  the  total  externally  applied  load  P which  will 

T 

cause  incipient  lateral  buckling.  The  only  loading  con- 
sidered is  distributed  load  in  the  negative  y direction. 

The  type  of  loading  (uniform,  triangular  or  whatever)  is 
specified  beforehand  by  a dimensionless  function  f(z)  such 
that 

w(z)  = (R/L3)  f(z)  (2) 

where  the  load  shape  function  w(z)  has  dimensions  of  force 
per  unit  length.  Concentrated  loads  and  concentrated  moments 
can  be  approximated  by  appropriate  large  local  variations  in 
f(z).  The  function  f(z)  will  be  denoted  by  x5  when  it 
appears  later.  The  constant  R has  dimensions  of  force  times 
distance  squared.  The  magnitude  of  R is  arbitrary  and  is 
chosen  at  the  convenience  of  the  user. 

The  load  function  has  been  defined.  We  are  interested 
in  determining  the  multiplying  factor  Q such  that  lateral 
buckling  does  not  take  place  if  the  actual  loading  is  less 
than  Qw(z)  and  does  take  place  if  the  actual  loading  exceeds 
Qw(z).  The  entire  purpose  of  the  program  developed  from 
this  analysis  is  to  determine  this  multiplier  Q.  In  order 
to  obtain  nondimensionalized  equations,  it  is  convenient  to 
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introduce  P the  total  applied  load  for  incipient  buckling, 
defined  by 

L 

P = Q J w(z)  dz  (3) 

0 

We  will  also  employ  the  normalization 

L 

f(z)  dz  = L (4) 

0 

so  that  Q = PL2/R  (5  ) 

For  the  special  case  where  the  load  contributions  in  the 
negative  y direction  are  canceled  by  the  load  contributions 
in  the  positive  y direction  the  following  alternate  form  of 
equation  3 is  applied. 

L 

P = J Q | w ( z ) | dz  (6) 

0 

The  load  may  be  applied  above,  on,  or  below  the  cen- 
troidal  axis  as  specified  by  the  load  eccentricity,  e(z), 
cf.  figure  6 (page  20). 

For  notational  and  programming  convenience  a series  of 
dimensionless  functions  x are  defined.  The  first  six  are 
used  to  specify  beam  properties  and  loading.  A seventh  is 
simply  a function  having  constant  unit  value. 
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xl 

s 

EI1/R 

(7) 

x2 

s 

EI2/R 

(8) 

x3 

s 

C/R 

(9) 

X4 

= 

RL2/C1 

(10) 

x5 

= 

w(z)  L3/R  = f ( z ) 

(11) 

x6 

= 

e ( z )/L 

(12) 

x7 

S 

1 

(13) 

The  nondimensionalizing  scheme  used  to  define  the  x 
functions  is  applied  throughout  the  problem  development. 

The  dimensionless  forms  for  vertical  and  horizontal  deflec- 
tion are 


and 


xlO  = v/QL 
6 = u/L 


(14) 

(15) 


where  8 is  a dimensionless  function  which  ultimately  will  be 
constructed  from  several  functions. 

We  will  use  the  dimensionless  parameter 


C = z/L 


(16  ) 


to  specify  axial  position.  Accordingly  we  will  think  of  the 

x functions  as  functions  of  £ , thus  x5  = x5 ( c ) • Differentia- 

d Y 

tion  with  respect  to  z and  £ are  indicated  as  (p-  = Y* ) and 
dY 

(jt-  = Y)  respectively.  The  two  operations  are  related  by 
the  equation 


d 1 d 

3T  “ r 3"c 


(17) 
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The  developments  which  follow  require  considerable  inte- 
gration. Two  notational  devices  are  used  to  facilitate  both 
the  readability  and  the  transcription  of  equations,  Integra 
tion  is  signified  by  p,  thus 

C 

px3  = J x3(Od£  • 

0 


Evaluation  of  any  function  at  the  right  end  of  the  beam,  i.e 
for  £ s l,  is  indicated  by  the  prefix  *.  Several  examples 
illustrating  this  notation  are 

1 

*pxl  = J xl(C)dC  (18) 

0 

*x5x9  = x5 ( 1 ) • x9 ( 1 ) (19) 


*px5x6pxll 


1 


*x2  px8  = x2(l) 

1 

= j x5 ( C)x6 ( C ) [ 
0 


x8(c)d? 


C 


J xll(0)d0]d? 
0 


(20) 


(21) 


(In  the  last  equation  0 is  a dummy  variable.) 
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Figure  3.  Diagram  of  forces  and  moments  in  the  vertical 
plane . 

Our  analysis  is  restricted  to  beams  which  are  uncon- 
strained except  at  one  or  both  ends.  There  are  six  possible 
constraints  at  each  end  each  of  which  is  modeled  as  a spring. 
The  spring  constant  k describing  such  a constraint  can 
actually  vary  from  zero  (perfectly  free,  no  constraint  at 
all)  to  infinity  (perfect  or  ideal  constraint)  inclusive . 
Therefore  to  avoid  mathematical  difficulties  the  analysis 
employs  dimensionless  spring  parameters,  designated  by  the 
symbol  a (with  an  appropriate  identifying  suffix).  For  a 
torsion  spring,  such  as  that  corresponding  to  the  moment 
klv'(O)  shown  in  figure  3 (page  15),  the  relations  between 
kl  and  al  are 

al  = kl/ (kl  + R/L)  (22) 

kl  = Ral/L(l  - al)  (23) 

Note  that  al  = 0 corresponds  to  kl  = 0 (perfect  freedom) 
and  al  = 1 corresponds  to  kl  = ® 


(perfect  constraint). 


Similarly,  for  a linear  spring  such  as  spring  number  3,  the 
relations  are 

a3  s k3/(k3  ♦ R/L3)  (24) 

k3  = Ra3/L3 ( 1 - a3)  (25) 

Figure  5 (page  19)  shows  that  the  point  of  application  of  a 
linear  spring  is  not  limited  to  the  centroidal  axis.  The 
vertical  deflection  constraints  may  be  applied  at  distances 
cl  and  c2  above  the  centroidal  axis,  and  the  horizontal  de- 
flection constraints  may  be  applied  at  distances  c3  and  c4 
above  the  centroidal  axis.  The  positive  sense  is  as  indi- 
cated in  figure  5,  and  the  dimensionless  forms  of  these  dis- 
tances are 


G1 

= cl/L 

(26) 

G2 

= c2/L 

(27) 

G3 

= c3/L 

(28) 

G4 

= c4/L 

(29) 

The  physical  and  geometrical  parameters  used  in  the  prob- 
lem have  now  been  defined.  They  will  be  used  in  the  struc- 
tural analysis  developed  in  the  next  section. 
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III.  STRUCTURAL  ANALYSIS 

There  are  three  structural  deformations  examined  in  the 
analysis.  They  are  vertical  bending,  horizontal  bending, 
and  torsion. 

The  vertical  bending  problem  is  disjoint  from  the  remain- 
der of  the  analysis,  and  it  is  solved  first.  The  notation 
used  in  defining  vertical  bending  will  not  be  shown  because 
the  problem  is  simple  and  familiar.  It  is  sufficient  to  say 
that  a set  of  four  simultaneous  nonhomogeneous  equations  are 
formed,  representing  the  constraint  conditions  at  each  end 
which  pertain  to  bending  in  the  vertical  (yz)  plane.  This 
system  is  then  solved  by  a standard  algebraic  analysis,  and 
the  result  is  used  to  create  three  dimensionless  functions 
for  shear,  moment,  and  deflection. 


Shear 

x8 

= V/P 

(30) 

Moment 

x9 

= M/PL 

(31) 

Deflection 

xlO 

= v/QL 

(32) 

The  remainder  of  the  analysis  involves  studying  the  com- 
bined horizontal  bending  and  torsion  problem.  Figure  4 
(page  18)  shows  the  positive  sense  for  the  constant  horizon- 
tal shear  H,  the  linearly  varying  horizontal  bending  moment 
N,  and  horizontal  displacement  u. 


J 


17 


k7Cu(0)+cl+(0)3  H kaCu<l>*e2*(l>D 


Figure  4.  Diagram  showing  forces  and  moments 
in  the  horizontal  plane. 

The  appropriate  dimensionless  variables  are 


8 

= u/L  ; 

u = 

L8 

(33) 

B4 

= H/P  ; 

H = 

PB4  = B4P 

(34) 

n 

= N/PL; 

N = 

PLn 

(35) 

Dimensionless  shear  B4  is  one  of  ultimately  eight  un- 
known scalar  quantities  (Bl,  B2  , . . . . , B8)  which  will  be  in- 
troduced in  the  process  of  determining  Q.  The  alternate 
forms  indicated  in  the  preceding  equations  are  intended  to 
familiarize  the  reader  with  notational  forms  which  will  be 
freely  employed  in  later  developments. 

The  illustration  of  the  rotated  cross  section  in  figure 
2 (page  10)  shows  that  the  moment  causing  bending  about  the 
n axis  is  N cos  <J>  - M sin  4> , which  reduces  to  N - M<f>  by  small 
angle  approximation.  From  elementary  bending  theory  we  have 

EI2  u"  = N - <J>M  (36) 
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and  by  substitution  from  equation  7 it  may  be  written  as 


B = Qx2 (n  - 0x9)  (37) 

Examination  of  the  simple  statics  of  figure  4 (page  18) 
yields 

Hz  + N(0)  = N(z)  (38) 

In  dimensionless  form  this  becomes 

B4C  + n(0)=  n ( C ) (39) 

« 

and  may  be  rewritten  as 

n = B4px7  + B5x7  (40) 

n(0)  = B5  (41) 


r 


— ■ 


Figures  5 and  6 (page  20)  show  the  forces  and  moments 
used  to  develop  the  torsion  moment  equation 


T(z)  = T( 0 )-H[ v( z )-v(0 ) ] + V(z)u(z)-V(0)u(0) 
z 

+ J w(z)  [u(z)+e(z)<J»(z)  ]dz  (42) 

0 


Introducing  the  dimensionless  function 

t = T/PL  (43) 


and  using  additional  dimensionless  functions  previously  de- 
fined, equation  31  may  be  written  as 

t = t(0)-QB4[xl0-xl0(0)]+x8  6-x8(0)8(0  )+px5  ( 0+x6<{> ) (44) 


V(z),down 


• T(0) 


T(U 


wfm^V4(n>  v(L>r^4(l  

I.  LEFT  END 


c.  RIGHT  END 


Figure  6.  Diagram  of  left,  central,  and  right 
portions  for  torsion  analysis. 
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Reference  (1)  develops  the  following  third  order  equa- 
tion which  accounts  for  the  warping  rigidity  of  the  beam. 


. C1$'M  = C<T  + T - u'M  (45) 

In  dimensionless  form  this  may  be  written  as 

’♦  = x4(Qt  - Qx9 6 + x3$)  (46) 

The  function  x4  has  previously  been  defined  as 

x4  = RL2/C1  (47) 

Cl  is  the  coefficient  of  the  highest  (i.e.,  third)  order 
derivative  in  equation  45,  and  if  Cl  becomes  small,  very 
troublesome  mathematical  difficulties  are  encountered,  cf. 
Reference  (3).  Accordingly,  the  analysis  for  the  case  Cl 
= 0 cannot  be  obtained  from  that  for  Cl  i 0 , given  in  refer- 
ence (1),  except  with  the  greatest  of  difficulty.  It  is 
for  this  reason  that  the  separate  analysis  represented  by 
this  thesis,  was  undertaken.  The  analysis  herein,  to  this 
point,  is,  as  has  been  said  before,  essentially  identical 
to  that  in  reference  (1).  From  this  point  on  there  are 
essential  differences. 

The  present  analysis  presumes  that  Cl  = 0 , so  that  equa- 
tion 35  takes  the  form 

l = Qx3_1 ( x9 6 - t)  (48) 

The  inverse  of  x3,  x3“\  occurs  frequently  in  what 
follows.  The  term  is  clumsy  in  this  form  and  will  be  re- 
designated as 
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(49) 


Therefore  equation  48  is  written  as 

♦ = QY(x9 8 - t)  (50) 

Two  of  the  end  constraint  conditions  of  the  original 
problem,  namely  those  relating  to  warping  constraints  at  the 
end,  are  not  available  for  and  do  not  appear  in  the  present 
specialized  problem. 

The  original  problem  had  six  constraints  at  each  end. 

The  present  problem  has  five  constraints  at  each  end.  Four 
constraint  conditions  were  employed  in  solving  the  problem 
of  bending  in  the  vertical  plane.  Thus  a total  of  six  con- 
straint conditions  remain  at  this  stage  of  the  development. 
The  first  four  result  from  the  constraints  against  horizon- 
tal deflection  and  end  rotations  in  the  x-z  plane.  The  last 
two  result  from  constraints  against  rotation  about  the  z 


axis . 

a5 8(0 ) - Q(l-a5)  n(0)  = 0 (51) 

a68(l)  + Q(l-a6)  n(l)  = 0 (52) 

a7[ 8(0 )+Gl$(0 ) ] + QB4(l-a7)  = 0 .(53) 

a8  [ 6 ( 1 ) +G2 4> ( 1 ) ] - QB4(l-a8)  = 0 (54) 

a9$(0)+Q(l-a9)[t(0)+G3x8(0)<f>(0)-B4Gl]  = 0 (55) 

al0$(l  )-Q(  1-alO  ) [t(l)+G4x8(l)<j>(l)-B4G2  ] = 0 (56) 


Although  this  notation  appears  formidable,  it  is  con- 
venient, explicit,  and  unequivocal.  The  terms  are 
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sufficiently  complicated  that  no  matter  what  notation  might 
be  used  careful  attention  to  detail  is  required. 

There  is  no  need  to  demonstrate  the  complete  development 
of  the  six  boundary  equations,  but  for  illustrative  purposes 
the  derivations  of  two  typical  constraint  equations,  51  and 
53,  will  be  given. 

The  constraint  against  rotation  of  the  left  end  about  a 


vertical  (y)  axis,  as  shown  in  figure  6 (page  20)  results 
in  the  following  equation. 

k5u' (0)  - NCO)  = 0 (57) 

Employing  the  dimensionless  variable  8 = u/L  and  n = N/PL, 
using  equation  17,  writing  the  spring  constant  k5  in  terms 
of  the  spring  parameter  a5 , viz, 

k5  = Ra5/L(l-a5)  (58) 

and  recalling  equation  5,  it  is  easy  to  obtain  the  equation 

a5B(0)-Q(l-a5)n(0)  = 0 (59) 

This  is  equation  51. 

Equation  53  defines  the  constraint  against  horizontal 
deflection  of  the  left  end  of  the  beam.  The  deflection  is 
u(0)  + cl$(0),  see  figure  5 (page  19),  where  the  second  term 
results  from  permitting  the  spring  to  be  attached  at  a dis- 
tance cl  above  the  centroidal  axis.  The  constraint  equation 
is 


where  as  previously  defined  H = PB4 , cl  = LG1,  and 

k7  = a7R/L3 (l-a7 ) (61) 

Using  the  dimensionless  quantities  previously  introduced 
it  is  a simple  matter  to  put  equation  60  in  the  form  ex- 
hibited as  equation  53. 

The  structural  analysis  has  now  been  completed  and  is  re 
presented  by  equations  37,  44,  50,  51,  52,  53,  54,  55,  and 

1 

56.  The  mathematical  method  of  solution  will  be  addressed 


IV.  MATHEMATICAL  ANALYSIS 


Before  going  into  the  details  of  the  analysis,  we  will 
try  to  describe  it  in  general  terms. 

An  expression  is  formed  for  $ which  involves  an  unknown 
function  (xll  in  what  follows)  the  initial  form  of  which  is 
arbitrarily  assumed  except  for  the  normalizations 

xll ( 0 ) = 0 (62) 

*p | xll | = 1 (63) 

The  corresponding  function  d>  is  obtained  by  integration  and 
equation  37  is  used  to  get  8.  This  is  followed  by  two  more 
integrations  to  get  8 and  8.  These  are  employed  in  the  tor- 
sion equation  44  to  get  the  function  t.  Finally  this  is 
used  in  equation  48  to  get  a new  expression  for  4> . Each  in- 
tegration introduces  a new  unknown  scalar  quantity.  In  one 
way  or  another  a total  of  eight  such  unknown  scalar  constants 
Bl,  B2,...,B8  are  introduced.  The  six  constraint  equations, 
plus  two  more  which  will  be  introduced  in  what  follows,  pro- 
vide a set  of  eight  linear  homogeneous  equations  in  these  un- 
known B's  which  thus  leads  to  an  eigenproblem  of  eighth 
order.  The  parameter  Q appears  in  the  coefficients  to  the 
first  and  second  powers.  The  quadratic  eighth  order  eigen- 
problem is  transformed  to  a linear  sixteenth  order  eigenprob- 
lem. Since,  generally,  both  zero  and  infinite  eigenvalues 
are  contained  in  the  solution,  a special,  relatively  new 
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algorithm,  the  Q2  method,  capable  of  dealing  with  such  sys- 
tems, is  employed  for  the  solution  of  this  eigenproblem. 

The  appropriate  eigenvalue  is  chosen  and  the  corresponding 
eigenvector  provides  the  solution  for  the  unknown  B's  so 
that  the  new  xll  can  be  uniquely  determined.  This  new  xll 
is  used  in  place  of  the  original  xll  and  the  process  is  re- 
peated until  there  is  convergence.  Mo  proof  is  offered 
that  indeed  convergence  must  ultimately  take  place;  because 
of  the  complexity  of  the  problem  su  .h  a proof  might  be  very 
difficult  to  establish.  However,  experience  with  the  pro- 
cedure described  here  indicates  that  convergence  may  indeed 
be  expected.  The  details  of  this  procedure  are  now  given  in 
what  follows. 

The  function  # is  initially  specified  as 

= 31x11  + 32x7  (64) 

where  xll  is  an  arbitrarily  assumed  function  satisfying  equa- 
tions 62  and  63  and  where  31  and  B2  are  the  third  and  fourth 
of  eight  unknown  constants  that  will  be  employed.  The  first 
was  B4  from  equation  34  and  the  second  was  35  from  equation 

41. 

Integrating  equation  64  yields 

= Blpxll  + 32px7  + B3x7  (65) 

and  when  equations  40  and  65  are  substituted  into  equation  3" 
we  arrive  at 
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6 = -BlQx2x9pxll  - B2Qx2x9px7  - B3Qx2x9 


+B4Qx2px7  + B5Qx2  (66) 

Integrating  twice  more,  we  get 

8 = -BlQpx2x9pxll  - B2Qpx2x9px7  - B3Qpx2x9 

+B4Qpx2px7  + B5Qpx2  + B6Qx7  (67) 

B = -BlQppx2x9pxll  - B2Qppx2x9px7  - B3Qppx2x9 

+B4Qppx2px7  + B5Qppx2  + B6Qpx7  + B7Qx7  (68) 

For  convenience  the  next  step  is  to  rewrite  equation  44 
as 

t = qB  + px5x64>  - QB4X  - Qx8(0)B7  + QB8 

where  t(0)  = QB8 

X = xlO  - xlO(O) 
q = px5  + x8 

The  symbol  q denotes  an  operator  acting  on  the  function  B, 
not  a function  which  is  multiplied  by  8. 

When  the  preceding  equations,  evaluated  at  C=0  or  C= 1 as 
the  case  may  be,  are  substituted  in  boundary  conditions  equa- 
tions 51  through  56,  a system  of  six  linear  simultaneous 
homogeneous  equations  in  eight  unknown  constants  B is  obtained. 

In  the  previous  section,  which  addressed  the  structural 
analysis,  equations  51  and  53  were  derived  for  illustrative 
purposes.  Now  for  illustrative  continuity  these  two  equa- 
tions will  be  used  to  show  the  method  of  formulation  of  the 
six  linear  simultaneous  homogeneous  equations. 


(69) 


(70) 

(71) 
(72  ) 
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Recalling  that  equation  51  is 

a5  6(0 ) - Q(l-a5)  n(0)  = 0 
and  substituting 

6(0)  = QB6 
n (0  ) = B5 

we  now  factor  out  Q leaving 

-B5(l-a5)  + B6a5  = 0 


(73) 

(74) 


(75) 


The  second  illustrative  equation,  53,  is 

a7  [ 6(0  ) + Gl<f>  ( 0 ) ] + QB4(l-a7)  = 0 


and  upon  substituting 

4>  (0 ) = B3  . (76  ) 

6(0)  = B7Q  (77) 

we  arri ve  at 

B3Gla7  + B4Q(l-a7)  + B7Qa7  = 0 (78) 

If  equations  67  and  69  are  substituted  into  equation  50 
the  following  "new"  <j>  is  created. 
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i = Y[B1(-Q2x9px2x9pxll  + Q‘ qppx2x9pxll  - Qpx5x6pxll) 

2 2 

+ B2(-Q  x9px2x9px7  + Q qppx2x9px7  - Qpx5x6px7) 

+ B3(-Q2x9px2x9  + Q2qppx2x9  - Qpx5x6) 

+ B4(Q2x9px2px7  + Q2X  - Q2qppx2px7) 

+ B5(Q2x9px2  - Q2qppx2) 

+ B6(Q2x9  - Q2qpx7) 

+ B8(-Q2)]  (79) 

(The  coefficient  of  B7  vanishes  identically.) 

To  complete  the  eigensystem  involving  the  eight  unknown 
B's,  two  additional  equations  are  needed. 

The  seventh,  equation,  consistent  with  equation  62,  is 

new  <j>(0)  - B2  = 0 (80) 

The  eighth  equation  is 

<*P*>0ld  - * 0 <91) 

and  is  arrived  at  by  observing  that  if  i given  in  equation  64 

and  created  from  an  assumed  xll,  were,  somehow,  a "correct" 

• • 

then  it  and  the  recovered  in  equation  79  would  be  iden- 
tically equal  to  each  other.  The  integration  in  equation  81 
achieves  a certain  "smoothing"  or  averaging  which  has  been 
found  to  be  preferable  to  simply  equating  the  new  and  old 
functions  at  an  arbitrary  value  of  C. 

The  coefficients  of  the  eigensystem  equations  generally 
involve  Q to  the  zeroth,  first,  and  second  powers  so  that 
the  system  may  be  exhibited  in  the  form 

(D  + EQ  + FQ2)B  = 0 (82) 
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The  nonzero  elements  of  matrices  D,  E,  and  F are 


D15  = a5-l,  D16  = a5 , D21  = -a6*px2x9pxll , 

D22  = -a6*px2x9px7 , D2 3 = -a6*px2x9,  D24  = a6*px2px7-a6  + l , 
D25  = a6*px2-a6+l,  D26  = a6 , D33  = Gla7 , 

D41  = a8G2*pxll,  D42  = a8G2 , D43  = a8G2 , 

D53  = a9 , D61  = alO*pxll,  D62  = alO , D63  = alO , 

D72  = -1,  D81  = -*pxll,  D82  = -1,  E34  = l-a7, 

E37  = a7,  E41  = -a8 Appx2x9pxll , 

E42  = -a8ftppx2x9px7 , E43  = -a8,vppx2x9, 

E44  = a8ftppx2px7+a8-l , E45  = a8*ppx2, 

E46  = a8 , E47  = a8,  E53  = ( l-a9 )G3x8 (0  ) , 

E54  = (a9-l)Gl,  E61  = (alO-1 ) [ ftpx5x6pxll+G4*x8 ( 1 )pxll ] , 

E62  = (alO-1) [*px5x6px^  + G4*x8(l)], 

E63  = (alO-1)  [ ,'px5x6+G4,ltx8  ( 1 ) ] , E64  = (l-alO)G2, 

E81  = -*pYpx5x6pxll,  E82  = -ftpYpx5x6px7 , 

E83  = -ftpYpx5x6 , F58  = l-a9 , 

F61  = ( 1-alO ) Aqppx2x9pxll , F62  = ( 1-alO ) *qppx2x9px7  , 

F63  = ( 1-alO ) rtqppx2x9 , F64  = ( 1-alO ) C *X-*qppx2px7 ) , 

F65  = (al0-l)*qppx2 , F66  = (alO-1 )*qpx7 , 

F68  = alO-1 , F76  = x9(0)Y(0),  F78  = -Y(0), 

F81  = -ftpYx9px2x9pxll  + *pYqppx2x9pxll , 

F82  = -ApYx9px2x9px7  + ,vpYqppx2x9px7 , 

F83  = -ftpYx9px2x9  + ApYqppx2x9, 

F84  = ftpYx9px2px7  + ftpYX  - ftpYqppx2px7 , 

F85  = ftpYx9px2  - 'vpYqppx2  , 

F86  = *PYx9  - *pYqpx7, 

F88  = -*pY 
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(For  reference  purposes  all  the  preceding  equations  bear  the 
number  84 . ] 

This  eighth  order  quadratic  eigensystem  may  be  expanded 
to  a sixteenth  order  linear  eigensystem  in  the  form 


VZ  = QWZ  (85) 

where  V = I (86) 

Li  oj 


and  I and  0 respectively  represent  the  eighth  order  unit  and 
null  matrices. 

The  eigenproblem  has  now  been  defined  in  a linear  form, 
but  because  the  matrices  V and  W may  be  singular  not  just 
any  method  will  satisfactorily  solve  the  problem.  We  have 
chosen  the  QZ  method  described  in  reference  (3).  From  the 
sixteen  eigenvalues  produced  the  smallest  positive  eigenvalue 
is  selected.  The  last  eight  elements  of  the  associated  eigen- 
vector are  the  desired  constants  B1  through  B8. 

Using  equation  79  a new  $ is  now  formed  and  from  it  a 
new  xll  is  created  by  means  of  the  following  normalization 
process . 


xll 


new 


[♦ 


new 


- $ (0)]/[*p|<fr 

new  ' 1 F 1 rnew 


- ♦ (0)  I ] 

new  1 1 


(89) 
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The  iteration  process  is  now  repeated  until  satisfactory 


I 


convergence  is  attained.  The  converged  smallest  positive 
eigenvalue  Q is  now  used  to  determine  the  buckling  load  from 
either  equation  3 or  equation  6. 


r 

i 

I 
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V.  COMPUTER  IMPLEMENTATION  OF  THEORY 

To  this  point  the  functions  x have  been  treated  as  con- 
tinuous functions  for  use  in  an  analytic  solution.  However, 
it  is  obvious  that  solving  the  problem,  even  for  simple 
functions,  by  analytic  means  is  prohibitively  difficult. 

To  overcome  this  problem  a numerical  method  of  solution  was 
employed.  The  beam  was  divided  into  a finite  number  of 
equal  sections.  The  functions  x were  then  redefined  as  vec- 
tors with  Fortran  designations  x(l,I),  x(2,I)...,  etc. 
where  I = 1,  2,  3,...,N.  The  elements  of  any  such  vector 
represent  function  values  at  the  N points  of  subdivision  of 

i , 

| 

the  beam  into  equal  subsections  of  length  L/CN-1).  The 
vector  x(10,I)  was  used  twice,  first  to  represent  xlO  and 
then  to  store  X of  equation  71. 

The  vectors  in  this  form  readily  lend  themselves  to 
manipulation.  The  simple  subroutines  developed  to  perform 
the  operations  are,  with  the  exception  of  the  integration 
scheme,  not  worthy  of  note  here. 

Trapezoidal  integration  was  first  used  during  the  debug- 

ging  of  the  raw,  untested  program.  This  method  gives  satis- 

2 

factory  results  with  an  error  of  0(h  ),  where  h is  the  in- 
terval length.  Once  it  was  established  that  the  program 
would  give  satisfactory  solutions  the  refinement  process  was 
begun.  A higher  order  numerical  integration  scheme  was  em- 
ployed next.  It  is  described  by  Milne  in  reference  (5), 
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and  results  in  an  error  0(h5). 

The  QZ  subroutine  is  a canned  program  in  the  library  of 
the  Naval  Postgraduate  School  Computer  Center.  It  was  trans- 
scribed  from  reference  (6),  with  very  little  change,  by 
Mr.  Roger  Hilleary.  To  date,  it  has  been  the  only  program 
that  has  given  satisfactory  solutions  to  the  present  eigen- 
value problem. 

The  complete  theory  described  in  the  preceding  sections 
of  this  thesis  has  been  implemented  in  the  form  of  a sub- 
routine LATBRO,  a complete  listing  of  which  is  contained  in 
Appendix  A. 

The  user  must  write  a main  program  which  supplies  dimen- 
sioning statements , number  of  intervals , specifications  for 
vectors  xl  through  x6 , specifications  for  the  constants  G, 
and  the  ten  spring  parameters.  The  following  is  an  example 
of  such  a main  program.  It  was  used  to  solve  what  is  listed 
as  Case  6,  Table  1,  page  37. 

IMPLICIT  REAL* 8 (A-H,0-Z) 

REAL* 8 X(30 ,101) ,ALFA(12 ) ,G ( 4 ) , S , ONE , ZERO 
INTEGER  KP(IO) 

COMMON  X ,ALFA,G ,N, KP 
N=101 

ZERO=0 . D+0 
0NE= 1. D+0 
DO  1 1 = 1, N 
X(1 , I )=0NE 
X(2 ,I)=0NE 
X ( 3 , 1 ) = 0NE 
X( 5 , 1 ) = ZERO 
X ( 6 , 1 ) = ZERO 

1 CONTINUE 
X(5 , 1 )=0NE 
DO  2 1=1,4 

2 G(I)=ZER0 
DO  3 1=1,12 


3 ALFA( I ) =ZERO 
ALFA C 2 )=ONE 
ALFA(4)=0NE 
ALFA (6 )=ONE 
ALFA( 8 )=ONE 
ALFA (10)= ONE 
KP ( 3 ) = 1 

CALL  LATBRO(S) 

STOP 

END 


, i 
' ; 


I 

f 
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VI.  TESTING  AND  VERIFICATION 

As  mentioned  in  the  previous  section,  the  numerical  in- 
tegration methods  that  have  been  employed  have  their  short- 
comings . 

If  the  integrand  is  "smooth"  either  the  trapezoidal  or 
the  higher  order  Milne  method  may  be  used.  The  latter  is 
preferable  because  of  its  greater  accuracy. 

Both  methods  introduce  errors,  however,  when  there  are 
steps,  impulses,  or  doublets  to  be  integrated.  We  have 
found  that  the  Milne  and  trapezoidal  methods  can  be  used  to 
solve  problems  with  a concentrated  load.  The  error  depends 
linearly  on  the  reciprocal  of  the  number  of  elements  into 
which  the  beam  is  divided  and  extrapolation  may  be  used  to 
obtain  an  excellent  result.  As  yet  we  do  not  have  experi- 
ence with  concentrated  moment  loadings  (doublets).  This 
subject  is  discussed  in  more  detail  in  reference  (1). 

The  following  table  shows  the  cases  that  have  been  tes- 
ted to  date.  In  each  of  these  cases  all  the  G's  and  the 
eccentricity  x6  were  zero.  Each  of  these  has  been  success- 
ful. 

The  values  for  Q shown  in  the  last  column  of  the  table 
were  obtained  with  x2  = 1,  x3  = 1,  and  by  simple  scaling 
analysis  lead  to  total  load 

P = Q /dfel2/L2  (89) 
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Cases  of  Uniform  Beams  Analyzed  by  Subroutine  LATBRO 
(x6  = 0,  G1  = G2  = G3  = G4  = 0) 


1 

j J 

for  lateral  buckling.  The  values  reported  agree  with  values 
to  be  found  in  reference  (2). 

The  following  are  details  of  several  cases. 

Case  1 

As  shown,  it  is  a cantilever  beam  fully  restrained  at  the 
right  end.  The  loading  is  unifoim. 


Figure  7.  Illustration  of  side  view  of  Casa  1 beam. 


The  subroutine  produced  a converged  value  of  Q of  12.854. 

Timoshenko's  result  is  12.85;  cf.  page  261  ref  (2). 

Case  2 

This  is  simply  Case  1 turned  end  for  end.  It  was  used 
to  test  additional  equations  and  also  resulted  in  a converged 
critical  load  output  of  12.854  which  is  identical  to  the  re- 
suit  for  Case  1,  as  it  should  be. 

!j 
1 1 
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Case  5 


The  converged  critical  load  output  was  16.936  and  Timo- 
shenko's result  is  16.94 j cf.  page  269  ref  (2). 


/\ 
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Figure  8.  Top  and  side  views  Case  5 beam. 

In  several  of  the  cases  listed  in  the  table  above,  the 
loading  does  not  interact  with  one  or  more  of  the  springs, 
so  that  the  corresponding  spring  parameters  could  have  an 
arbitrary  value  without  changing  the  results. 
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The  subroutine  LATBRO  in  its  present  form  should  be  able 
to  solve  beam  problems  with  beams  that  have  uniform  or  vari- 
able cross  sections,  with  loading  on  the  centroidal  axis, 
and  with  spring  constraints  attached  at  the  centroidal  axis. 
Loading  that  involves  concentrated  moments , as  mentioned  be- 
fore, has  not  yet  been  tested. 

The  program  was  designed  to  solve  problems  where  loading 
and  spring  constraints  may  act  away  from  the  centroidal  axis 
It  is  recommended  that  follow  on  testing  begin  with 

1.  loading  off  the  centroidal  axis. 

2.  spring  constraints  attached  away  from  the  centroidal 
axis . 

3.  concentrated  moment  loading. 

It  is  expected,  with  the  exception  of  minor  programming 
errors,  that  problems  with  any  combination  of  the  above  will 
be  readily  solved  by  the  subroutine. 
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APPENDIX  A 


Listing  of  LATBRO  and  Adjunct  Subroutines 
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